Setup

source("helpers.R")

dir.create("results/DEA", showWarnings = F, recursive = T)

knitr::opts_chunk$set(fig.width = 10, dpi = 300, results = "hold", fig.show = "hold")
# Heatmaps

hm_cluster_rows <- TRUE # Genes
hm_cluster_cols <- TRUE # Samples
hm_scale_by_row <- TRUE
hm_max_rows <- 20

heatmap.colors <- colorRamp2(c(-2, 0, 2), c("darkblue", "white", "darkred"))

# Volcano plot

vp_lfc_limit <- 5

Undifferentiated

Supernatant

data <- read_DIA_report("data/MS1MS2/20240701_082307_ASC_SN_TLR10_Light_Report.tsv")

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")  
plot_missing(data_MS1.df, data_MS2.df)

keep <- filter_too_many_missing(data_MS1.df, data_MS2.df, full.groups = 1)

data.filtered <- data[rownames(data) %in% keep,]


data_MS1.filtered.df <- assay(data.filtered, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.filtered.df <- assay(data.filtered, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.filtered.df, data_MS2.filtered.df)

data.log2 <- data.filtered
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

meanSdPlot(assay(data.log2, "MS1"))
meanSdPlot(assay(data.log2, "MS2"))

data.norm <- data.filtered
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

MNAR_MS1 <- get_MNAR(data_MS1.norm.df)
MNAR_MS2 <- get_MNAR(data_MS2.norm.df)
MNAR <- MNAR_MS1 | MNAR_MS2

plot_missing_heatmap(data_MS1.norm.df, data_MS2.norm.df, MNAR, colData(data.norm),
                     colors_columns = list(Celltype=c("ASC_WT" = "lightblue",
                                                      "ASC_TLR10_LOV" = "lightgreen"),
                                           Condition=c("Light" = "orange",
                                                       "Dark" = "midnightblue")),
                     colors_rows = list(MNAR=c("TRUE"="darkgreen",
                                               "FALSE"="darkmagenta")))

plot_missing_density(data_MS1.norm.df, data_MS2.norm.df)

data.imputed <- data.norm

assays(data.imputed) %<>% lapply(MsCoreUtils::impute_mixed, !MNAR, "MLE", "MinDet")
## Loading required namespace: norm
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
data_MS1.imputed.df <- assay(data.imputed, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.imputed.df <- assay(data.imputed, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
## Iterations of EM: 
## 1...2...3...4...5...6...
## Iterations of EM: 
## 1...2...3...4...5...
plot_pca(data.imputed, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDarkInWT = "ASC_WT.Light - ASC_WT.Dark", 
  LightVsDarkInTLR10 = "ASC_TLR10_LOV.Light - ASC_TLR10_LOV.Dark",
  TLR10VsWTInDark = "ASC_TLR10_LOV.Dark - ASC_WT.Dark",
  TLR10VsWTInLight = "ASC_TLR10_LOV.Light - ASC_WT.Light",
  Diff = "(ASC_TLR10_LOV.Light - ASC_TLR10_LOV.Dark) - (ASC_WT.Light - ASC_WT.Dark)"
)

fit <- fit_DEqMS_model(data.imputed, contrast_list)

VarianceBoxplot(fit)

ASC Wildtype Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "WT")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_WT_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_WT_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC Wildtype Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 189 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TLR10LOV Light vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "TLR10_LOV")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_TLR10LOV_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_TLR10LOV_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TLR10LOV Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 262 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TLR10LOV vs. Wildtype Dark

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_TLR10LOV_vs_WT_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_TLR10LOV_vs_WT_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TLR10LOV vs. Wildtype Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 192 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TLR10LOV vs. Wildtype Light

current_contrast <- 4
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Light")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_TLR10LOV_vs_WT_Light_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_TLR10LOV_vs_WT_Light_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows,
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TLR10LOV vs. Wildtype Light Supernatant", vp_lfc_limit)
## Warning: ggrepel: 64 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC Interaction (TLR10LOV Light vs. Dark) vs. (Wildtype Light vs. Dark)

current_contrast <- 5
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, ]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_Interaction_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_Interaction_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, cluster_rows = hm_cluster_rows,
             cluster_cols = hm_cluster_cols, scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC Interaction Supernatant", vp_lfc_limit)

Whole Cell Lysate

data <- read_DIA_report("data/MS1MS2/20240701_082255_ASC_WCL_TLR10_Light_Report.tsv")

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.df, data_MS2.df)

keep <- filter_too_many_missing(data_MS1.df, data_MS2.df)

data.filtered <- data[rownames(data) %in% keep,]


data_MS1.filtered.df <- assay(data.filtered, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.filtered.df <- assay(data.filtered, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.filtered.df, data_MS2.filtered.df)

data.log2 <- data.filtered
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

meanSdPlot(assay(data.log2, "MS1"))
meanSdPlot(assay(data.log2, "MS2"))

data.norm <- data.filtered
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

MNAR_MS1 <- get_MNAR(data_MS1.norm.df)
MNAR_MS2 <- get_MNAR(data_MS2.norm.df)
MNAR <- MNAR_MS1 | MNAR_MS2

plot_missing_heatmap(data_MS1.norm.df, data_MS2.norm.df, MNAR, colData(data.norm),
                     colors_columns = list(Celltype=c("ASC_WT" = "lightblue",
                                                      "ASC_TLR10_LOV" = "lightgreen"),
                                           Condition=c("Light" = "orange",
                                                       "Dark" = "midnightblue")),
                     colors_rows = list(MNAR=c("TRUE"="darkgreen",
                                               "FALSE"="darkmagenta")))

plot_missing_density(data_MS1.norm.df, data_MS2.norm.df)

data.imputed <- data.norm

assays(data.imputed) %<>% lapply(MsCoreUtils::impute_mixed, !MNAR, "MLE", "MinDet")
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
data_MS1.imputed.df <- assay(data.imputed, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.imputed.df <- assay(data.imputed, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
## Iterations of EM: 
## 1...2...3...4...5...
## Iterations of EM: 
## 1...2...3...4...5...
plot_pca(data.imputed, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDarkInWT = "ASC_WT.Light - ASC_WT.Dark", 
  LightVsDarkInTLR10 = "ASC_TLR10_LOV.Light - ASC_TLR10_LOV.Dark",
  TLR10VsWTInDark = "ASC_TLR10_LOV.Dark - ASC_WT.Dark",
  TLR10VsWTInLight = "ASC_TLR10_LOV.Light - ASC_WT.Light",
  Diff = "(ASC_TLR10_LOV.Light - ASC_TLR10_LOV.Dark) - (ASC_WT.Light - ASC_WT.Dark)"
)

fit <- fit_DEqMS_model(data.imputed, contrast_list)

VarianceBoxplot(fit)

ASC Wildtype Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "WT")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()

print(paste("Found", nrow(data.sig.df), "differentially expressed proteins."))

write.csv(res, file = "results/DEA/ASC_WT_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_WT_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
## [1] "Found 162 differentially expressed proteins."
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC Wildtype Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 121 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TLR10LOV Light vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "TLR10_LOV")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_TLR10LOV_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_TLR10LOV_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TLR10LOV Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 36 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TLR10LOV vs. Wildtype Dark

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_TLR10LOV_vs_WT_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_TLR10LOV_vs_WT_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TLR10LOV vs. Wildtype Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 462 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC TLR10LOV vs. Wildtype Light

current_contrast <- 4
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Light")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_TLR10LOV_vs_WT_Light_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_TLR10LOV_vs_WT_Light_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows,
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC TLR10LOV vs. Wildtype Light Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 121 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ASC Interaction TLR10LOV Light vs. Dark vs. Wildtype Light vs. Dark

current_contrast <- 5
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, ]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ASC_Interaction_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ASC_Interaction_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ASC_WT" = "lightblue",
                                                       "ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ASC Interaction Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 91 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC

Supernatant

data <- read_DIA_report("data/MS1MS2/20240701_082311_ADIPO_ASC_SN_TLR10_Light_Report.tsv")

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.df, data_MS2.df)

keep <- filter_too_many_missing(data_MS1.df, data_MS2.df)

data.filtered <- data[rownames(data) %in% keep,]


data_MS1.filtered.df <- assay(data.filtered, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.filtered.df <- assay(data.filtered, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.filtered.df, data_MS2.filtered.df)

data.log2 <- data.filtered
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")


plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

data.norm <- data.filtered
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

MNAR_MS1 <- get_MNAR(data_MS1.norm.df)
MNAR_MS2 <- get_MNAR(data_MS2.norm.df)
MNAR <- MNAR_MS1 | MNAR_MS2

plot_missing_heatmap(data_MS1.norm.df, data_MS2.norm.df, MNAR, colData(data.norm),
                     colors_columns = list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                      "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                           Condition=c("Light" = "orange",
                                                       "Dark" = "midnightblue")),
                     colors_rows = list(MNAR=c("TRUE"="darkgreen",
                                               "FALSE"="darkmagenta")))

plot_missing_density(data_MS1.norm.df, data_MS2.norm.df)

data.imputed <- data.norm

assays(data.imputed) %<>% lapply(MsCoreUtils::impute_mixed, !MNAR, "MLE", "MinDet")
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
data_MS1.imputed.df <- assay(data.imputed, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.imputed.df <- assay(data.imputed, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
## Iterations of EM: 
## 1...2...3...4...5...6...
## Iterations of EM: 
## 1...2...3...4...5...6...7...
plot_pca(data.imputed, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDarkInWT = "ADIPO_ASC_WT.Light - ADIPO_ASC_WT.Dark", 
  LightVsDarkInTLR10 = "ADIPO_ASC_TLR10_LOV.Light - ADIPO_ASC_TLR10_LOV.Dark",
  TLR10VsWTInDark = "ADIPO_ASC_TLR10_LOV.Dark - ADIPO_ASC_WT.Dark",
  TLR10VsWTInLight = "ADIPO_ASC_TLR10_LOV.Light - ADIPO_ASC_WT.Light",
  Diff = "(ADIPO_ASC_TLR10_LOV.Light - ADIPO_ASC_TLR10_LOV.Dark) - (ADIPO_ASC_WT.Light - ADIPO_ASC_WT.Dark)"
)

fit <- fit_DEqMS_model(data.imputed, contrast_list)

VarianceBoxplot(fit)

ADIPO ASC Wildtype Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "WT")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_WT_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_WT_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC Wildtype Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 130 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC TLR10LOV Light vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "TLR10_LOV")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_TLR10LOV_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_TLR10LOV_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC TLR10LOV Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 467 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC TLR10LOV vs. Wildtype Dark

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC TLR10LOV vs. Wildtype Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 585 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC TLR10LOV vs. Wildtype Light

current_contrast <- 4
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Light")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Light_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Light_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows,
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC TLR10LOV vs. Wildtype Light Supernatant", vp_lfc_limit)
## Warning: ggrepel: 401 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC Interaction (TLR10LOV Light vs. Dark) vs. (Wildtype Light vs. Dark)

current_contrast <- 5
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, ]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_Interaction_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_Interaction_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, cluster_rows = hm_cluster_rows,
             cluster_cols = hm_cluster_cols, scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC Interaction Supernatant", vp_lfc_limit)

Whole Cell Lysate

data <- read_DIA_report("data/MS1MS2/20240626_081714_ADIPO_ASC_WCL_TLR10_Light_Report.tsv")

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.df, data_MS2.df)

keep <- filter_too_many_missing(data_MS1.df, data_MS2.df)

data.filtered <- data[rownames(data) %in% keep,]


data_MS1.filtered.df <- assay(data.filtered, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.filtered.df <- assay(data.filtered, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.filtered.df, data_MS2.filtered.df)

data.log2 <- data.filtered
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")


plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

data.norm <- data.filtered
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

MNAR_MS1 <- get_MNAR(data_MS1.norm.df)
MNAR_MS2 <- get_MNAR(data_MS2.norm.df)
MNAR <- MNAR_MS1 | MNAR_MS2

plot_missing_heatmap(data_MS1.norm.df, data_MS2.norm.df, MNAR, colData(data.norm),
                     colors_columns = list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                      "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                           Condition=c("Light" = "orange",
                                                       "Dark" = "midnightblue")),
                     colors_rows = list(MNAR=c("TRUE"="darkgreen",
                                               "FALSE"="darkmagenta")))

plot_missing_density(data_MS1.norm.df, data_MS2.norm.df)

data.imputed <- data.norm

assays(data.imputed) %<>% lapply(MsCoreUtils::impute_mixed, !MNAR, "MLE", "MinDet")
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
data_MS1.imputed.df <- assay(data.imputed, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.imputed.df <- assay(data.imputed, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
## Iterations of EM: 
## 1...2...3...4...
## Iterations of EM: 
## 1...2...3...4...5...
plot_pca(data.imputed, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDarkInWT = "ADIPO_ASC_WT.Light - ADIPO_ASC_WT.Dark", 
  LightVsDarkInTLR10 = "ADIPO_ASC_TLR10_LOV.Light - ADIPO_ASC_TLR10_LOV.Dark",
  TLR10VsWTInDark = "ADIPO_ASC_TLR10_LOV.Dark - ADIPO_ASC_WT.Dark",
  TLR10VsWTInLight = "ADIPO_ASC_TLR10_LOV.Light - ADIPO_ASC_WT.Light",
  Diff = "(ADIPO_ASC_TLR10_LOV.Light - ADIPO_ASC_TLR10_LOV.Dark) - (ADIPO_ASC_WT.Light - ADIPO_ASC_WT.Dark)"
)

fit <- fit_DEqMS_model(data.imputed, contrast_list)

VarianceBoxplot(fit)

ADIPO ASC Wildtype Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "WT")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_WT_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_WT_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC Wildtype Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 47 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC TLR10LOV Light vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "TLR10_LOV")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_TLR10LOV_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_TLR10LOV_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC TLR10LOV Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 66 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC TLR10LOV vs. Wildtype Dark

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC TLR10LOV vs. Wildtype Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 121 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC TLR10LOV vs. Wildtype Light

current_contrast <- 4
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Light")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Light_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_TLR10LOV_vs_WT_Light_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows,
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC TLR10LOV vs. Wildtype Light Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 90 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

ADIPO ASC Interaction (TLR10LOV Light vs. Dark) vs. (Wildtype Light vs. Dark)

current_contrast <- 5
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, ]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/ADIPO_ASC_Interaction_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/ADIPO_ASC_Interaction_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("ADIPO_ASC_WT" = "lightblue",
                                                       "ADIPO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "ADIPO ASC Interaction Whole Cell Lysate", vp_lfc_limit)

OSTEO ASC

Supernatant

data <- read_DIA_report("data/MS1MS2/20240701_082302_OSTEO_ASC_SN_TLR10_Report.tsv")

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.df, data_MS2.df)

keep <- filter_too_many_missing(data_MS1.df, data_MS2.df)

data.filtered <- data[rownames(data) %in% keep,]


data_MS1.filtered.df <- assay(data.filtered, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.filtered.df <- assay(data.filtered, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.filtered.df, data_MS2.filtered.df)

data.log2 <- data.filtered
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")


plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

data.norm <- data.filtered
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

MNAR_MS1 <- get_MNAR(data_MS1.norm.df)
MNAR_MS2 <- get_MNAR(data_MS2.norm.df)
MNAR <- MNAR_MS1 | MNAR_MS2

plot_missing_heatmap(data_MS1.norm.df, data_MS2.norm.df, MNAR, colData(data.norm),
                     colors_columns = list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                      "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                           Condition=c("Light" = "orange",
                                                       "Dark" = "midnightblue")),
                     colors_rows = list(MNAR=c("TRUE"="darkgreen",
                                               "FALSE"="darkmagenta")))

plot_missing_density(data_MS1.norm.df, data_MS2.norm.df)

data.imputed <- data.norm

assays(data.imputed) %<>% lapply(MsCoreUtils::impute_mixed, !MNAR, "MLE", "MinDet")
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
data_MS1.imputed.df <- assay(data.imputed, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.imputed.df <- assay(data.imputed, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
## Iterations of EM: 
## 1...2...3...4...5...6...
## Iterations of EM: 
## 1...2...3...4...5...6...
plot_pca(data.imputed, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDarkInWT = "OSTEO_ASC_WT.Light - OSTEO_ASC_WT.Dark", 
  LightVsDarkInTLR10 = "OSTEO_ASC_TLR10_LOV.Light - OSTEO_ASC_TLR10_LOV.Dark",
  TLR10VsWTInDark = "OSTEO_ASC_TLR10_LOV.Dark - OSTEO_ASC_WT.Dark",
  TLR10VsWTInLight = "OSTEO_ASC_TLR10_LOV.Light - OSTEO_ASC_WT.Light",
  Diff = "(OSTEO_ASC_TLR10_LOV.Light - OSTEO_ASC_TLR10_LOV.Dark) - (OSTEO_ASC_WT.Light - OSTEO_ASC_WT.Dark)"
)

fit <- fit_DEqMS_model(data.imputed, contrast_list)

VarianceBoxplot(fit)

OSTEO ASC Wildtype Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "WT")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_WT_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_WT_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC Wildtype Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 43 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC TLR10LOV Light vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "TLR10_LOV")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_TLR10LOV_Light_vs_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_TLR10LOV_Light_vs_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC TLR10LOV Light vs. Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 18 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC TLR10LOV vs. Wildtype Dark

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Dark_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Dark_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC TLR10LOV vs. Wildtype Dark Supernatant", vp_lfc_limit)
## Warning: ggrepel: 392 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC TLR10LOV vs. Wildtype Light

current_contrast <- 4
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Light")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Light_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Light_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows,
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC TLR10LOV vs. Wildtype Light Supernatant", vp_lfc_limit)
## Warning: ggrepel: 345 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC Interaction (TLR10LOV Light vs. Dark) vs. (Wildtype Light vs. Dark)

current_contrast <- 5
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, ]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_Interaction_SN_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_Interaction_SN_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, cluster_rows = hm_cluster_rows,
             cluster_cols = hm_cluster_cols, scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC Interaction Supernatant", vp_lfc_limit)

Whole Cell Lysate

data <- read_DIA_report("data/MS1MS2/20240701_082250_OSTEO_ASC_WCL_TLR10_Report.tsv")

data_MS1.df <- assay(data, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
data_MS2.df <- assay(data, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.df, data_MS2.df)

keep <- filter_too_many_missing(data_MS1.df, data_MS2.df)

data.filtered <- data[rownames(data) %in% keep,]


data_MS1.filtered.df <- assay(data.filtered, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.filtered.df <- assay(data.filtered, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
plot_missing(data_MS1.filtered.df, data_MS2.filtered.df)

data.log2 <- data.filtered
assays(data.log2) %<>% lapply(log2)

data_MS1.log2.df <- assay(data.log2, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.log2.df <- assay(data.log2, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")


plot_intensity_boxplot(data_MS1.log2.df, data_MS2.log2.df)

data.norm <- data.filtered
assays(data.norm) %<>% lapply(normalize_matrix, "vsn")

data_MS1.norm.df <- assay(data.norm, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.norm.df <- assay(data.norm, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

plot_intensity_boxplot(data_MS1.norm.df, data_MS2.norm.df)

meanSdPlot(assay(data.norm, "MS1"))
meanSdPlot(assay(data.norm, "MS2"))

MNAR_MS1 <- get_MNAR(data_MS1.norm.df)
MNAR_MS2 <- get_MNAR(data_MS2.norm.df)
MNAR <- MNAR_MS1 | MNAR_MS2

plot_missing_heatmap(data_MS1.norm.df, data_MS2.norm.df, MNAR, colData(data.norm),
                     colors_columns = list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                      "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                           Condition=c("Light" = "orange",
                                                       "Dark" = "midnightblue")),
                     colors_rows = list(MNAR=c("TRUE"="darkgreen",
                                               "FALSE"="darkmagenta")))

plot_missing_density(data_MS1.norm.df, data_MS2.norm.df)

data.imputed <- data.norm

assays(data.imputed) %<>% lapply(MsCoreUtils::impute_mixed, !MNAR, "MLE", "MinDet")
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
## Imputing along margin 1 (features/rows).
data_MS1.imputed.df <- assay(data.imputed, "MS1") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")

data_MS2.imputed.df <- assay(data.imputed, "MS2") %>%
  as.data.frame() %>%
  rownames_to_column("Gene.Name")
## Iterations of EM: 
## 1...2...3...4...5...
## Iterations of EM: 
## 1...2...3...4...5...6...
plot_pca(data.imputed, scale = T, plot_all = T, maxPC = 4)

contrast_list <- c(
  LightVsDarkInWT = "OSTEO_ASC_WT.Light - OSTEO_ASC_WT.Dark", 
  LightVsDarkInTLR10 = "OSTEO_ASC_TLR10_LOV.Light - OSTEO_ASC_TLR10_LOV.Dark",
  TLR10VsWTInDark = "OSTEO_ASC_TLR10_LOV.Dark - OSTEO_ASC_WT.Dark",
  TLR10VsWTInLight = "OSTEO_ASC_TLR10_LOV.Light - OSTEO_ASC_WT.Light",
  Diff = "(OSTEO_ASC_TLR10_LOV.Light - OSTEO_ASC_TLR10_LOV.Dark) - (OSTEO_ASC_WT.Light - OSTEO_ASC_WT.Dark)"
)

fit <- fit_DEqMS_model(data.imputed, contrast_list)

VarianceBoxplot(fit)

OSTEO ASC Wildtype Light vs. Dark

current_contrast <- 1
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "WT")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_WT_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_WT_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC Wildtype Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 7 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC TLR10LOV Light vs. Dark

current_contrast <- 2
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "TLR10_LOV")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_TLR10LOV_Light_vs_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_TLR10LOV_Light_vs_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC TLR10LOV Light vs. Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 66 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC TLR10LOV vs. Wildtype Dark

current_contrast <- 3
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Dark")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Dark_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Dark_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC TLR10LOV vs. Wildtype Dark Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 252 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC TLR10LOV vs. Wildtype Light

current_contrast <- 4
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, 
                         str_detect(colnames(data.imputed), "Light")]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Light_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_TLR10LOV_vs_WT_Light_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows,
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC TLR10LOV vs. Wildtype Light Whole Cell Lysate", vp_lfc_limit)
## Warning: ggrepel: 237 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

OSTEO ASC Interaction (TLR10LOV Light vs. Dark) vs. (Wildtype Light vs. Dark)

current_contrast <- 5
res <- outputResult(fit, coef_col = current_contrast) %>%
  as_tibble() %>%
  mutate(qvalue = qvalue(P.Value)$qvalues,
         sca.qvalue = qvalue(sca.P.Value)$qvalues) %>%
  select(!c(adj.P.Val, sca.adj.pval))

res.sig <- res %>%
  filter(sca.qvalue < 0.05 & abs(logFC) >= 0.58)

qvals.sig <- res.sig$sca.qvalue
names(qvals.sig) <- res.sig$gene

data.sig <- data.imputed[rownames(data.imputed) %in% res.sig$gene, ]

data.sig.df <- assay(data.sig, "MS2") %>% as.data.frame()
  

write.csv(res, file = "results/DEA/OSTEO_ASC_Interaction_WCL_unfiltered.csv",
          row.names = FALSE)

write.csv(res.sig, file = "results/DEA/OSTEO_ASC_Interaction_WCL_filtered.csv",
          row.names = FALSE)
ggplot(data = res, aes(x = sca.P.Value)) + 
  geom_histogram(binwidth = 0.025)

plot_heatmap(data.sig, heatmap.colors, list(Celltype=c("OSTEO_ASC_WT" = "lightblue",
                                                       "OSTEO_ASC_TLR10_LOV" = "lightgreen"),
                                            Condition=c("Light" = "orange",
                                                        "Dark" = "midnightblue")),
             qvalues = qvals.sig, title = "Rel. LogIntensity", max_rows = hm_max_rows, 
             cluster_rows = hm_cluster_rows, cluster_cols = hm_cluster_cols, 
             scale_by_row = hm_scale_by_row)

plot_volcano(res, "OSTEO ASC Interaction Whole Cell Lysate", vp_lfc_limit)